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ABSTRACT 

A cryogenic achromatic half-wave plate (HWP) for submillimetre astronomical 
polarimetry has been designed, manufactured, tested, and deployed in the Balloon- 
borne Large- Aperture Submillimeter Telescope for Polarimetry (BLASTPol). The de- 
sign is based on the five-slab Pancharatnam recipe and it works in the wavelength 
range 200-600 /im, making it the most achromatic HWP built to date at submillime- 
tre wavelengths. The frequency behaviour of the HWP has been fully characterised at 
room and cryogenic temperatures with incoherent radiation from a polarising Fourier 
transform spectrometer. We develop a novel empirical model, complementary to the 
physical and analytical ones available in the literature, that allows us to recover the 
HWP Mueller matrix and phase shift as a function of frequency and extrapolated to 
4K. We show that most of the HWP non-idealities can be modelled by quantifying one 
wavelength-dependent parameter, the position of the HWP equivalent axes, which is 
then readily implemented in a map-making algorithm. We derive this parameter for a 
range of spectral signatures of input astronomical sources relevant to BLASTPol, and 
provide a benchmark example of how our method can yield unprecedented accuracy 
on measurements of the polarisation angle on the sky at submillimetre wavelengths. 

Key words: instrumentation: polarimeters — techniques: polarimetric — balloons 
— magnetic fields — polarisation — submillimetre. 



1 INTRODUCTION 

Star formation occurs in the cold interstellar medium (ISM) 
where the gas is mostly found in molecular form. A small 
fraction, typically 10 -6 , of gas particles ionised by cosmic 
rays provides strong coupling between the cold gas and the 



ambient interstellar magnetic field. Thus magnetic fields 
might play a crucial role in the evolution of star-forming 
clouds, perhaps controlling the rate at which stars form 
and eve n determining the masses of stars (see review ar- 
ticles bv lCrutche3l2004lMaKee fc OstrikeifeOOTt) . Although 
the imp ortance of magneti c fields was suggested as early as 
1956 by iMe stel fc Spitzen . it has not yet been possible to 
clearly establish their influence on Giant Molecular Clouds 
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(GMCs) and star formation. This lack of understanding 
is caused by the difficultly in observing Galactic magnetic 
fields on the spatial scales relevant to the star-forming pro- 
cesses, especially within obscured molecu lar clouds (see e.g., 
ICrutcher et al]|2004l ; IWhittet et aTll200Sl ). 

The only direct method for determining the strength of 
the magnetic field is to measure the splitting of mo lecular 
lines through the Zeeman effect (e.g.. Ic"rutcherl [l999). This 
has been carried out successfully for a number of molecular 
cloud cores, though the t echnique is difficult an d is limited to 
the brightest of regions jCrutcher et al.lll999l ). In addition, 
this method only yields the component of the magnetic field 
along the line of sight, and so can only be used to determine 
the average field strength for a statistical sample of cores. 

One promising alternative method for probing Galac- 
tic magnetic fields is to observe clouds with a far- 
infrared /submillimetre (FIR/submm) polarimeter (e.g., 
iHildebrand et al ] [2003; IWard-Thomoson et all Exx]). By 
tracing the linearly polarised thermal emission from dust 
grains aligned with respect to the local magnetic fields, we 
can estimate the direction of the plane-of -the-sky compo - 
nent of the field within the cloud (see e.g.. lLazarian|[20CM) . 
and its strength via the Chandrasekhar & Fermi (CF; Il953l ) 
technique, provided that ancillary measurements of the tur- 
bulent motion velocity are available. 

Ground-based obse r vations with the SCU BA polarime- 
ter (jMurrav et al.lll997l : iGreaves et al1l2003l ) and the Sub- 
millimeter Pol_arhneto forAntarctic Remote Observations 
(SPARO; iNovak et ail 120031 ) show that the submm emis- 
sion from prestellar cores and GMCs is indeed polarised 
to a few perc ent dWard- Thompson et al ] l2000l : iLi et all 
2006). Planck (|Ade et al.ll201ll ) will provide coarse resolu- 
tion (FWHM ~5') submm polarimetry maps of the entire 
Galaxy. The Atacama Large Millimeter/ submillimeter Ar- 
ray (ALMA; IWootten fc Thompson! 120091 ) will provide sub- 
arcsecond resolution millimetre (mm) and submm polarime- 
try, capable of resolving fields within cores and circumstellar 
disks, but will not be sensitive to cloud-scale fields. 

The Balloon-borne Large- Apertu re Submillimete r Tele- 
scope for Polarimetry (BLASTPol; iMarsden et alj 120081 ; 
iFissel et al.l 120101 : iPascale et al.ll2012l ). with its arcminute 
resolution, is the first submm polarimeter to map the large- 
scale magnetic fields within molecular clouds with high sen- 
sitivity and mapping speed, and sufficient angular resolu- 
tion to observe into the dense cores (~0.1pc). BLAST- 
Pol mapped the polarised dust emission over a wide range 
of column densities corresponding to Ay > 4 mag, yielding 
hundreds to thousands of independent polarisation pseudo- 
vectors per cloud, for a dozen clouds. Moreover, the polari- 
metric observations of BLASTPol at 250, 35 0, and 500 nm 
complement those plan ned for SCUBA-2 (|Bastien et alj 
l2005l ; lHolland et al.l200^ ) at 850 /im, as BLASTPol's smaller 
pixel count is almost completely compensated by the in- 
creasing flux density at shorter wavelengths. In particular, 
BLASTPol has better sensitivity to degree-scale polarised 
emission. Core maps to be obtained using SCUBA-2 can be 
combined with those produced by BLASTPol to trace mag- 
netic structures in the cold ISM from scales of 0.01 pc out to 
5 pc, thus providing a much needed bridge between the large- 
area but coarse-resolution polarimetry provided by Planck 
and the high-resolution but limited field-of-view maps of 
ALMA. 



In January 2011, BLASTPol successfully completed its 
first 9.5-day flight over Antarctica. Ten science targets were 
mapped with unprecedented combined mapping speed and 
resolution. BLASTPol is scheduled to fly again in Decem- 
ber 2012 with the same apparatus to observe more targets 
and improve the statistics on GMCs and dark clouds. These 
observations will constitute an exciting dataset for studying 
the role played by magnetic fields in star formation. 

The BLASTPol linear polarisation modulation scheme 
comprises a stepped cryogenic achromatic half-wave plate 
(HWP) and photolithographed polarising grids placed in 
front of the detector arrays, acting as analysers. The grids 
are patterned to alternate the polarisation angle sampled by 
90° from bolometer-to-bolometer along the scan direction. 
BLASTPol scans so that a source on the sky passes along a 
row of detectors, and thus the time required to measure one 
Stokes parameter (either Q or U) is just equal to the separa- 
tion between bolometers divided by the scan speed. During 
operations, we carry out spatial scans at four HWP rotation 
angles spanning 90° (22.5° steps), allowing us to measure 
the other Stokes parameter through polarisation rotation. 

The use of a continuously rotating or stepped HWP 
as a polarisation modulator is a w idespread technique at 



mm and submm w a velengths (e.g.. iRenbarger et al 
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2009 1; [Johnson et al.ll2007l : ILi et al.ll2008l; IMatsumura et~aH 
I2OO9I : bryan et all l2010bF iDowell et all |2010| ). A thorough 
account of the HWP non-idealities and its inherent polari- 
sation systematics, especially for very achromatic designs, 
has become necessary as the accuracy and sensitivity of 
mm/submm instruments have soared in recent years. 

The literature offers numerous efforts to address, 
through simulations, the impact of the inevitable instru- 
mental systematic errors due to the polarisation modu- 
lation strategy in the unbiased recovery of the Stokes 
parameters Q, U on the sky, especially for Cosmic Mi- 
crowave Background (CMB) polarisation exp eriments (e.g., 
lO'Dea et all l2007l . l201ll : iBrown et aD 120091 ) . In addition, 
physical and analytical models have been developed to 
retrieve the frequency-dependent modulation function of 
achromatic HWPs and estimate the corrections d ue to 
non-flat source spectral indices (|Savini et all 120061 . 120091 : 
IMatsumura et al.ll2009l ). 

Nevertheless, little work has been published on in- 
corporating the measured HWP non-idealities in a data- 
analysi s pipeline and ultim ately in a map-making algo- 
rithm. iBrvan et all lj2010bT l derive an analytic model that 
parametrises the frequency-dependent non-idealities of a 
monochromatic HWP and present a map-making algo rithm 
that convincingly accounts for these. iBao et aD l|2012l ) care- 
fully simulate the impact of the spectral dependence of the 
polarisation modulation induced by an achromatic HWP on 
measurements of the CMB polarisation in the presence of 
astrophysical foregrounds, such as Galactic dust. However, 
both these works assume the nominal design values for the 
build parameters of the HWP plus anti-reflection coating 
(ARC) assembly. 

While this assumption is a reasonable one when no spec- 
tral measurements of the HWP as-built are available, sev- 
eral studies clearly show that the complex multi-slab crystal 
HWP and its typically multi-layer ARC are practically im- 
possible to manufacture exactly to the desired specifications. 
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In pa rticular, ISavini et all (|2006l . 120091 ) and iPisano et al.l 
(2006) caution against the finite precision to which the mul- 
tiple crystal substrates composing an achromatic H WP can 
be aligned relative to each other in the Pa ncharatnam (|l955l ) 
scheme. In addition, IZhang et all (2009) show how some of 
the design parameters in the ARC can slightly change during 
the bonding of la yers, achieved via a hot-pressing technique 
jAde et aljlSoOrj ). 

This work describes a novel empirical method that al- 
lows the reconstruction of the Mueller matrb|3 of a generic 
HWP as a function of frequency through spectral transmis- 
sion measurements of the HWP rotated by different angles 
with respect to the input polarised light. Not only does this 
method give complete and quantitative information on the 
measured spectral performance of the HWP, but it also pro- 
vides a direct avenue to accounting for the non-idealities of 
the HWP as-built in a map-making algorithm. This empiri- 
cal approach is applied to the BLASTPol HWP and will ulti- 
mately yield unprecedented accuracy on astronomical mea- 
surements of polarisation angles at submm wavelengths. 

The layout of this paper is as follows. In Section [2] we 
give an overview of the manufacturing process for BLAST- 
Pol's five-slab sapphire HWP. Section describes the spec- 
tral measurements, while Section [4] presents the empirical 
model as well as the main results of the paper. Finally, in 
Section [SI we describe the algorithm for the naive-binning 
map-making technique implemented by BLASTPol, which 
naturally accounts for the measured HWP non-idealities. 
Section [5] contains our conclusions. 



2 THE BLASTPol HALF- WAVE PLATE 

Wave plates (or retarders) are optical elements used to 
change the polarisation state of an incident wave, by induc- 
ing a predetermined phase difference between two perpen- 
dicular polarisation components. A (monochromatic) wave 
plate can be simply obtained with a single slab of uniax- 
ial birefringent crystal of specific thickness, which depends 
upon the wavelength and the index of refraction of the crys- 
tal. A birefringent crystal is characterised by four parame- 
ters, n e , n , Qf c , a a , the real part of the indices of refraction 
and the absorption coefficient (in cm - ) for the extraordi- 
nary and ordinary axes of the crystal. At a specific wave- 
length Ao, the phase shift induced by a slab is determined 
uniquely by its thickness d, and reads: 



. ... 2nd, 
/\ip (Ao ) = — — (n e 
Ao 



(1) 



Given the operating wavelength Ao, the required phase shift 
for the wave plate is achieved by tuning the thickness d. 

While monochromatic wave plates have been (and 
are still being) used in mm and submm a s tronomical po- 
larimeters (see e.g., iRenbarger et all 12004 iLi et all I200I 
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pendence of the phase shift on wavelength, expressed in 



L)..v /ell et al.ll201of i, the inherent de- 



Equation [T] constitutes an intrinsic limit in designing a po- 
larisation modulator that operates in a broad spectral range 
(i.e., is achromatic). 



2.1 Achromatic half-wave plate design 

Achromaticity is necessary for wave plates that are de- 
signed for use with mu lti-band bolometric receivers, such 
as BLASTPol PILOT dBernard et al]|2007D . or SCUBA-2 
(Bast ierTet al.| [2005: Savi ni et al.ll2009l ). To achieve a broad- 
band performance, multiple-slab solutions have been con- 
ceive d in the past (|Pancharatnam| [l955; Ti tle fc Rosenberel 
Il98ll ) to compensate and keep the phase shift approximately 
constant across the bandwidth, by stacking an odd number 
(usually 3 or 5) of birefringent substrates of the same ma- 
terial, which are rotated with respect to each other about 
their optica^] axes by a frequency-dependent set of angles. 

Achromatic wave plates have been designed and 
built for astronomical polarimeters at mm and submm 
wavelengths by many authors in the last decade 
llHananv et~al] 120051; IPisano et al. l l2006l ; ISavini et all 120061 . 
120091 ; iMatsumura et al.ll2009l ). following the Poincare spher e 
(PS) method first introduced by Pancharatnam l|l955l ). 
which we briefly recall here for completeness. The polari- 
sation state of a monochromatic wave in a given reference 
frame can be represented by a set of coordinates, latitude 
and longitude, on the PS that quantify, respectively, the el- 
lipticity angle (sin 2\ oc sin Aip) and the orientation angle 
of its major axis (tan2i/» oc cosAi^>). A linearly polarised 
state appears only on the equator (with ±Q and ±[7 at the 
four antipodes), while the left and right circularly polarised 
states (±V) lie at the north and south poles, respectively. 

Propagation of a wave through a single birefringent slab 
will rotate its polarisation state on the PS by an amount 
dependent on the relation between wavelength and thick- 
ness (Equation m, about an axis whose orientation depends 
upon the position of the optic axis of the wave plate with 
respect to the reference frame of the incoming polarisation 
state. Specifically, an ideal monochromatic half-wave plate 
produces one PS rotation of 180°, changing a linear polari- 
sation state to another one on the equator. 

When a polychromatic wave packet enters a multiple- 
slab HWP, the input polarised states of all wavelengths over- 
lap in a single point on the PS. After the rotation due to 
the first slab, the polarisation states of different wavelengths 
will be scattered along an arc on the PS, with separations 
that depend on the bandwidth AA of the wave packet. As 
anticipated, this effect can be compensated for by stacking 
together an odd number of birefringent slabs, rotated with 
respect to each other by a symmetric pattern of angles (a, /3, 
7, /3, an d a for 5 slabs) abou t th eir optical axes (as derive d 
by e.g., iPancharatnaml fl955l and iTitle fc Rosenberg! Il98ll ) . 
The various polarisation states regroup in a small area of the 
PS surface, thus achieving a nearly frequency-independent 
output polarisation state, within a certain AA. 



We adopt the Stokes formalism to represent the time- 

averaged polarisation state of electromagnetic ra diation; f or a r e- 
view of polarisation basics we refer the reader to lCollettl jl993h . 



2 We distinguish between "optic" axis of a crystal, i.e. the di- 
rection in which a ray of transmitted light experiences no bire- 
fringence, and "optical" axis, i.e. the imaginary line along which 
there is some degree of rotational symmetry in the optical system. 
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Figure 1. Exploded view of the BLASTPol HWP. We also show 
the two-layer anti-reflection coating described in Section 12.21 

We note that, strictly speaking, all the four parame- 
ters that characterise a crystal, n c , n , cv c , a , depend upon 
wavelength (as we will illustrate in detail for sapphire); in 
particular, the different frequency-dependence of the ordi- 
nary and extraordinary refraction indices enters Equation [1] 
in a non-trivial way, thus rendering the design of an achro- 
matic HWP increasingly difficult as AA broadens. 

Using the above PS method, we have designed a HWP 
for the BLASTPol instrument, which requires an extended 
frequency range to cover three adjacent 30% wi de sp ectral 
bands at 250, 350, and 500 /im. A Pancharatnam (1955J) five- 
slab design is chosen with axis orientations of a — 0°, /3 = 
26°, 7 = 90.3°, j3 = 26°, and a = 0°; these angles are 
op timised using t he ph ysical and analytical model developed 
by ISavini et al.l (|2006t ) for an achromatic HWP, whic h in 
turn is based on the work of iKennaugh fc Adachil (|196Gl ). In 
Fig. [T] we show an exploded view of the BLASTPol HWP 
assembly; to our knowledge, with its ~100% bandwidth, this 
is the most achromatic half-wave plate manufactured to date 
at submm wavelength^]. 

2.2 HWP manufacture 

In addition to the broad spectral range of operation, the 
BLASTPol HWP i s required to func tion at cryogenic tem- 
peratures (4K; see lFissel et aljfeoiol ) for two main reasons: 
(1) reduce the thermal emission from a warm element placed 
in the optical path, which would constitute a significant 
background load on the bolometers; and (2) reduce the losses 
in transmission due to absorption from the stack of five crys- 
tal substrates, which drops dramatically with temperature. 
The absorption in a crystal at FIR wavelengths is the result 
of the interactive coupling between the incident radiation 
and phonons - the thermally induced vibrations of the con- 
stituent atoms of the substrate crystal lattice. Because the 

3 The H WP designed for the mm - wave E and B Experiment 
(EBE X; iMatsumura et alj |2009| ; iReichborn-Kiennerud et al.l 
|201CT) is nominally slightly more achromatic, with a ~110% band- 
width. However, a comprehensive spectral characterisation of the 
final (as-flown) AR,-coated HWP assembly has yet to be pub- 
lished. 



phonon population is much reduced at low temperatures, 
cooling the crystal effectively reduces the absorption. 

The two obvious candidates (uniaxial birefringent) crys- 
tals are sapphire and quartz, because of their favorable op- 
tical properties in the FIR/submm. Sapphire is chosen over 
quartz due to its larger difference between ordinary and ex- 
traordinary refraction i ndices (An e -o ~ 0.32 fo r sapphire, 
and ~ 0.032 for quartz; [Loewenstem et aT, 1973), which im- 
plies a smaller thickness for the substrates (see Equation [1]). 
Since quartz and sapphire have a comparable level of ab- 
sorption at cry ogenic temperatures in th e wavelength range 
of 200-600 jiim (jLoewenstein et al1ll97St ). thinner substrates 
are desirable to minimise absorption losses. Nonetheless, the 
thin sapphire substrates chosen for the BLASTPol HWP do 
indeed show appreciable absorption, especially at the short- 
est wavelengths (250 /im band; see Section [3]) . 

The five slabs of the Pancharatnam (1955) design all 
have the same thickness. To cover the broad wavelength 
range of 200-600 /im, a substrate thickness is chosen to pro- 
duce a HWP at the central wavelength of the central band, 
350 /jm. By using the valu es of the refractive indice s for 
cold sapphire pub li shed by lLoewenstein et ail (|l973l ) and 
ICook fc Perkowitj l| 19851 . Anfl o Mm » 0.32), and imposing 
the required phase shift of 180° between the two orthogonal 
polarisations travelling through the plate, Equation Q] yields 
for the thickness of a single substrate a value ~0.547mm. 
The nearest available thickness on the market is 0.5 mm. A 
deviation of ~0.047mm from the desired thickness trans- 
lates to a departure of ^15° from the ideal phase shift of 
180° at 350 £tm, which is approximately what we measure 
(see Fig. 1221) . We briefly discuss the implications of this sys- 
tematic effect at the end of Section f4. 31 

The orientation of the optic axis on each sapphire sub- 
strate is determined with a polarising Fourier transform 
spectrometer (pFTS hereafter) , which is briefly described in 
Section [3.11 Each substrate is rotated between two aligned 
polarisers at the pFTS output until a maximum signal is 
achieved. The use of two polarisers avoids any complication 
from a partially polarised detecting system and any cross 
polarisation incurred from the pFTS output mirrors. The 
HWP is assembled by marking the side of each substrate 
with its reference optic axis and rotating each element ac- 
cording to the Pancharatnam design described in the pre- 
vious section. The stack of five carefully-oriented sapphire 
substrates, interspersed with 6 /im layers of polyethylene, are 
bonded together with a hot-pressing technique used in stan- 
dard FIR/submm filter production (lAde et al.ll2006h . The 
polyethylene has negligible effects on the final optical per- 
formance of the HWP, because when heated it seeps into the 
roughened surfaces of the adjacent substrates. 

In order to improve the robustness of the bond, the 
individual substrates are sandblasted with aluminium ox- 
ide (AI2O3) prior to fusion; this procedure dramatically im- 
proves the grip of the polyethylene between adjacent crystal 
surfaces. Careful cleansing and degreasing of all the crys- 
tal surfaces is required after sandblasting; in particular, we 
found trichloroethylene to be most effective in removing the 
traces of oily substances due to the sandblasting process. 

The thickness of the resulting stack (uncoated HWP) 
is 2.55 ± 0.01mm; its diameter is 100.0 ± 0.1mm. A two- 
layer broadband anti-reflection coating (ARC), necessary 
to maximise the in-band transmission of the HWP, is also 
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hot-pressed to the front and back surfaces of the assembled 
plate, again using 6 //m layers of polyethylene. The layer ad- 
jacent to the sapphire is an artificial dielectric metamaterial 
(ADM) com posed of metal-mes h patterned onto polypropy- 
lene sheets (|Zhang et al.l [2009h . while the outer layer is a 
thin film of polytctrafiuoroethylene (PTFE). The thickness 
of the final stack (coated HWP) is 2.80±0.01 mm. The diam- 
eter of the ARC is set to 88.0 ±0.f mm, slightly smaller than 
that of the HWP to avoid any contact betwe en the coating 
and the HWP mount (see iFissel et alj|20ich ; the ARC is 
bonded concentrically to the HWP and thus its diameter 
defines the optically-active area of the HWP. 

Because of the thermal expansion mismatch between 
the sapphire and the polypropylene, the HWP assembly has 
been cryogenically cycled numerous times prior to the flight 
to test the robustness of the bond at liquid helium tem- 
peratures. The HWP has been successfully installed in the 
BLASTPol cryogenic receiver and has flown from a balloon 
platform for about ten days, without delamination of the 
ARC or damage to the assembly. 

However, for cryogenic crystal HWPs larger in diam- 
eter than the BLASTPol one, the application of a metal- 
mesh ADM as an ARC has proven extrem ely challeng- 
ing. T herefore, extending previous work by iPisano et al.l 
(2008), we have recently designed and realised a pro- 
totype polypropylene-embedded met al-mesh broadban d 
achromatic HWP for mm wavelengths (|Zhang et al.ll201lT ); 
this will allow next generation experiments with large- 
aperture detector arrays to be equipped with large-format 
(>20cm in diameter) HWPs for broadband polarisation 
modulation. 



3 SPECTRAL CHARACTERISATION 
3.1 Introduction 

We fully characterise the spectral performance of the 
BLAS TPol HWP by using a pFTS of the Martin-Puplett 
119701 ) type. The source is an incoherent mercury arc lamp 
with an aperture of 10 mm, whose emission is well approx- 
imated by a blackbody spectrum at T e g ~ 2000 K; a low- 
pass filter blocks radiation from the source at wavelengths 
shorter than ~3.4 /im. The interferometer is equipped with a 
P10 beam divider, a P2 input polariser (at the source), and 
a P10 output polariser. The pFTS has a (horizontally) po- 
larised output focused beam with //3.5 or, in other words, 
a converging beam with angles 9 < 8°. 

As we will show in the next sections, the pFTS allows 
us to measure the HWP performance as a function of fre- 
quency and incoming polarisation state. Furthermore, be- 
cause of the strong dependence of the sapphire absorption 
coefficient on temperature, we measure the spectral response 
of the HWP both at room (Section 13. 2\ and cryogenic tem- 
peratures (~120K; Section r3.3|) . Ultimately, we want to re- 
trieve the frequency-dependent HWP Mueller matrix and 
phase shift, which, in turn, determine its spectral response 
and modulation efficiency. 
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Figure 2. Schematic drawing of the room-temperature spectral 
measurements setup. The horizontally polarised output of a pFTS 
feeds into a polyethylene lens that creates a quasi-parallel beam 
and is then refocused onto the horn aperture of the bolometric 
detector. Two polarisers alternatively parallel and perpendicular 
create the necessary polarisation selection for the "co-pol" and 
"cross-pol" sets of measurements. The arrows for PP1 and PP2 
indicate the selected polarisation, so that t he wire grid orient a- 
tions are perpendicular to the arrows, (from IZhang et af1l2011r t . 



3.2 Room-temperature measurements 

The schematic drawing of the room-temperature measure- 
ment configuration is shown in Fig. [2] In the following, we 
describe each element in sequential order from the polarised 
pFTS output to the detector system. 

In order to measure the HWP performance at near- 
normal incidence, we use a planar convex polyethylene lens 
(with focus at the position of the output pFTS image) to 
generate a quasi-parallel beam section; a second lens refo- 
cuses the beam onto the horn aperture of the detector sys- 
tem. The maximum range of incident angles is thus limited 
by the input source aperture (10 mm), a beam spread of only 
1.6°. This ensures an even illumination of the entire HWP 
optically-active MI Cel. clS i f it were inside the BLAS TPol op- 
tics box (see lFissel et al.ll2010l ; |Pascale et al.ll2012T ). 

The HWP is placed centrally in the collimated beam 
section between two P10 polarisers (PP1 and PP2), which 
are tilted by 45° with respect to the optical axis to avoid 
standing waves between the optical elements. This tilt in- 
troduces four ports that are optically terminated with a 
close to ideal blackbody, Eccosorb AN72 absorbed- The effi- 
ciency of these polarisers is separately determined to exceed 
99.8% over the range of frequencies of interest, with a cross- 
polarisation of less than 0.1%. The polarisers are initially 
aligned with respect to each other, with the grid wires ver- 
tical (thus selecting horizontal polarisation) with respect to 
the optical bench, in order to avoid any projection effect 
when tilted. 

Following the convention depicted in Fig. [2] measure- 
ments with aligned polarisers are referred to as "co-pol" 
transmission, T cp . As shown in the next section, the HWP 
has a complementary response when the output polariser 
PP2 (analyser) is rotated by 90° about the optical axis of 



4 We denote with Pn [/im] the period of a photolithographcd 
wire grid polariser, which has n/2 copper strips with n/2 gaps on 
a 1.5 /an mylar substrate. 



5 Emerson and Cuming, Microwave Products, http://www ec- 
cosorb. com/. 
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the system (i.e., horizontal wires, selecting vertical polarisa- 
tion); data taken with this configuration are also necessary 
to completely characterise the HWP, and are referred to as 
"cross-pol" transmission^, T xp . 

Common to both the warm and cold measurements is 
the requirement to position and rotate the HWP accurately 
with respect to its optical axis. When at room temperature, 
the HWP is held and rotated by a motorised rotating mount 
positioned centrally between the two tilted polarisers. The 
mount has a fixed orientation with respect to the optical 
axis of the system; we position it so that the collimated 
beam has normal incidence on the HWP (within 1°), and 
evenly illuminates its surface. The electronically-controlled 
rotating mount can rotate the HWP about its optical axis to 
obtain the polarisation modulation; the resolution of the dig- 
ital angular encoder on the rotation angle is 0.001° (a CAD 
drawing of the optical bench setup, incl uding the motorise d 
rotating mount, can be found in Fig. 1 of lPisano et al . 2006). 

Finally, the detecting system used is a 4.2 K liquid he- 
lium cooled indium antimonide (InSb) detector, which is 
cryogenically filtered to minimise photon noise. The spectral 
coverage of the data is thus defined by the cut-off frequency 
of the light collector waveguide (5 cm -1 ) and by two low- 
pass filters in the cryostat housing the bolometric detector 
(60 cm" 1 ) . We pay particular attention to ensure the absorp- 
tion of any diffracted or reflected stray radiation. Besides 
terminating all unused optical ports as described above, ad- 
ditional Eccosorb AN72 covers all the exposed metallic sur- 
faces close to the optical path. 

The rapid scan system records interferograms with a 
8 jj,m sampling interval over a 10 cm optical path difference, 
at a scan speed of lcms -1 ; this results in a Nyquist fre- 
quency of 625 cm -1 and a spectral resolution of 0.05 cm -1 . 

A first dataset is obtained in co-pol configuration by 
scanning the spectrometer in the absence of the HWP, which 
we refer to as the background spectrum. This dataset de- 
fines the pFTS reference spectral envelope, and it is the set 
against which all the following spectra are divided in order 
to account for the spectral features of the source, pFTS op- 
tics, detector system, and laboratory environment (i.e., wa- 
ter vapour). Subsequently, the HWP is inserted in between 
the tilted polarisers in co-pol configuration, and spectra are 
acquired at many different HWP rotation angles (resulting 
in a data cube). To enhance the spectral signal-to-noise ra- 
tio, each dataset at a given angle is obtained by computing 
the Fourier transform of an (apodised and phase-corrected) 
average of 60 interferograms with the mirror scanned in both 
the forward and backward directions. As anticipated, the re- 
sulting spectra are divided by the background dataset, which 
in turn is the average of three spectra, to obtain the trans- 
mission of the coated HWP alone as a function of frequency. 

Because these data are collected over several hours, the 
amount of water vapour in the room is likely to slightly 
change with time; we account for this by taking background 
spectra approximately every hour and dividing the HWP 
spectra taken within that hour only by the corresponding 
background dataset. Nevertheless, discernible residuals from 



6 We note that this definition of cross-pol may differ from other 
conve ntions adopted in the literature (e.g., that of iMasi et al.l 
120061. who operate without a HWP). 
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Figure 3. Synthetic transmission spectrum from an atmospheric 
model, in arbitrary units. Provided by Ade (2009, private com- 
munication). 



atmospheric features can still be appreciated in the final 
HWP spectra, especially at wavenumbers k > 30 cm -1 (the 
BLASTPol 250 /im band) . We use a synthetic atmospheric 
transmission spectrum (provided by Ade 2009, private com- 
munication; shown in Fig. [3} to correct the original spectra 
by concurrently scaling the amplitude of the most promi- 
nent features, which are due to water vapour. We find that 
while some of the spectra do not need any correction at 
all, others need to be corrected by as much as ^15%; the 
corrected spectra are shown in Fig. [4] where each line is a 
spectrum at a different rotation angle of the HWP (in the 
range 6 = 0°-332°). 

An ideal HWP modulates the polarisation at 4 8, there- 
fore in a complete revolution there are four maxima (and 
minima), two for each of the birefringent axes. The zero an- 
gle in this case coincides with the HWP maximum, which is 
the HWP angle at which we measure maximum total power 
on the detector; this of course includes signal outside of 
the HWP bands (in the range 5-60 cm" 1 ). As we show in 
Section 14.21 the position of the equivalent axes of the sap- 
phire plate stack (and hence the position of the HWP max- 
ima/minima) depends upon the wavelength. Therefore the 
HWP maxima (and minima) we assign while taking spec- 
tra are just rough approximations. Although we increase 
the angle sampling rate in the vicinities of a maximum or 
minimum, in order to fully characterise the HWP it is not 
necessary to take spectra exactly at its maxima or minima. 

Due to polarisation symmetry, no appreciable change 
should be observed in pairs of datasets taken at angles that 
are 180° apart. We verify that the experimental setup is 
symmetric with respect to the HWP rotation by comparing 
spectra taken, for instance, near the two maxima, at O^ 1 *** = 
[0°, 180°] and at O 1 ™* = [88°, 268°]. The fact that the curves 
are superimposed confirms that there are no artifacts arising 
from misalignments in the optical setup. 

Although we do correct for the residual contaminations 
due to atmospheric features, which mainly affect the shorter 
BLASTPol wavelengths, we cannot rule out the possibility 
that some of the spectral fringes may still be altered. Fur- 
thermore, and more importantly, these spectra show signif- 
icant in-band transmission loss due to the absorption from 
sapphire at room temperature, which becomes more promi- 



BLASTPol HWP 7 




17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 
Wavenumber [cm" 1 ] 

Figure 4. Measured co-pol transmission spectra of the coated 
BLASTPol HWP at room temperature. Each line is obtained at 
a different HWP rotation angle by computing the Fourier trans- 
form of an (apodised and phase-corrected) average of 60 inter- 
ferograms. The resulting spectra are corrected for residual con- 
taminations due to atmospheric features by using the synthetic 
spectrum shown in Fig. \3\ The solid vertical black lines show the 
approximate extent of the three BLASTPol bands. 



nent with increasing frequency. Because of these two reasons, 
we do not take cross-pol spectra at room temperature and 
repeat our measurements with the HWP in a vacuum cav- 
ity, at temperatures as low as currently possible with the 
experimental apparatus at our disposal. 



3.3 Cold measurements 

The experimental setup for spectral measurements of the 
cold HWP is substantially different than that described in 
the previous section, except for the radiation source and the 
main pFTS module. 

We position the HWP in a liquid nitrogen-cooled remov- 
able module retrofitted in the vacuum cavity at the output 
port of the pFTS; a photograph and a brief description of 
the module, which we refer to as "cold finger", are given 
in Fig. [5] After the roughly two hours needed for the cold 
finger module to thermalise, its base plate reaches temper- 
atures close to 77 K, while the HWP holder thermalises at 
about 120 K, despite the careful thermal insulation and link 
to the base plate. Other cryogenic tests conducted by bond- 
ing a thermometer at the center of a single slab of sapphire 
ensure that the temperature measured at the edge of an 
aluminium or copper holder closely matches that of the sap- 
phire substrate at its center. 

While maintaining a constant level of liquid nitrogen in 
the cold finger, we can characterise the spectral response of 
the cold HWP, by rotating it inside the vacuum cavity with 
a rotation angle resolution of 0.06°. In this configuration, 
the P10 output polariser of the pFTS acts as PP1 in the 
room temperature setup (see Fig. [2}, while a second P10 
polariser (analyser, acting as PP2) is installed at the exit 
port of the vacuum cavity. On the outside of the cavity, the 
cryostat housing the bolometric detector is connected with 
no air gaps to the exit port. This time we use a composite 
bolometer cooled to 1.5 K by pumping on the liquid helium 




Figure 5. Photograph of the "cold finger" module, which fits in 
the vacuum cavity at the output port of the pFTS. The central 
cylinder is hollow and has to be regularly replenished with liquid 
nitrogen to maintain the temperature of the HWP at ~120K. 
Aluminium insulation and a thick copper strap improve the ther- 
mal performance of the module. Two thermometers monitor the 
temperature at the bottom of the cylinder (base plate) at the 
edge of the copper HWP holder. The rotator is manually driven 
via a gear train and a vacuum-seal shaft leading to a manual 
knob outside the module. The resolution of the analogue encoder 
on the rotation angle is 0.06°. The presence of a thermometer on 
the rotating element prevents rotations greater than ~180°. 

bath; this detector is again cryogenically low-pass filtered at 
60 cm -1 to minimise photon noise. 

Over two days of measurements, we acquire data cubes 
for co-pol (Figs. [6] and [8} and cross-pol (Figs. [7]and|9]) trans- 
missions using exactly the same parameters as quoted in 
the previous section, except for the scan speed, which we 
increase to 2cms~ 1 to speed up the measurement process 
at no expense to the quality of the spectra. The background 
dataset is obtained in co-pol configuration by scanning the 
spectrometer in the absence of the cold finger. Because of the 
controlled environment in the vacuum cavity, our measure- 
ments are now much less susceptible to the external environ- 
ment; however, we repeat background scans at the very end 
of our measurement session to monitor drifts in the bolome- 
ter responsivity and other potential systematic effects. Prior 
to inserting the cold finger in the cavity, we characterise the 
instrumental cross-pol of this setup by rotating PP2 by 90° 
in cross-pol configuration and acquiring three spectra. By 
averaging these cross-pol spectra and dividing by the co-pol 
background, we measure a cross-pol level of 0.2% or less 
across the entire spectral range of interest (5-60 cm -1 ); wc 
include the resulting cross-pol spectrum in Figs. [6] and [7] 
(dark pink line). 

In the surfaces depicted in Figs. [8] and [9] slices of the 
data cube along the wavenumber axis constitute the mea- 
sured spectra for different angles of the HWP, while slices 
along the angle axis represent the modulation function of the 
wave plate at a given frequency or, more precisely, within a 
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Figure 6. Measured co-pol transmission spectra equivalent to 
those shown in Fig. g] but with the HWP at ~120K. 
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Figure 7. Measured spectra of the HWP at ~120K equivalent 
to those shown in Fig. [6] but for cross-pol transmission. 





Figure 9. Equivalent data cube to that shown in Fig. [8] but for 
cross-pol transmission. Note how the two surfaces are comple- 
mentarily in counterphase to each other. Each measured spectra 
(as shown in Fig. (TJl is a slice of the surface perpendicular to the 
angle axis. 

narrow band of frequencies denned by a combination of spec- 
tral resolution and the spectrometer's instrument response 
function. 

The features visible in all spectra (including those 
shown previously in Fig. Q are spectral fringes due to stand- 
ing waves generated inside the stack of dielectric substrates 
(even with a quasi-perfect impedance matching coating on 
the outer surfaces) ; the presence of several interspersed lay- 
ers of polyethylene enhances the amplitude of the fringes by 
introducing small amounts of absorption at every internal 
reflection. 

3.4 Uncertainties on the measured spectra 

Because we average 60 interferograms to obtain the final 
spectrum at each HWP position, the statistical uncertainty 
associated with the average on a single dataset is found to 
be negligible, as expected. Rather, we decide to average to- 
gether all the available background interferograms that are 
collected over one day of measurements, and take their sta- 
tistical dispersion as our estimate of the uncertainty associ- 
ated with all the spectra collected on that day. Because the 
thermodynamic conditions in the cavity under vacuum are 
not susceptible to changes in the external environment, this 
procedure allows us to account for drifts in the bolometer re- 
sponsivity and other potential systematic effects. We show in 
Fig. [10] the mean background spectra and the associated er- 
ror for the co-pol and cross-pol measurement sessions. These 
errors are used in SectionfJJto estimate the uncertainties on 
the HWP Mueller matrix coefficients. 



Figure 8. Data cube represented by a surface obtained by stack- 
ing a set of spectral co-pol transmissions of the HWP at different 
angles. Each measured spectra (as shown in Fig. |6} is a slice of 
the surface perpendicular to the angle axis. 



3.5 Modulation function and efficiency 

We can reduce the dependence on frequency of our data 
cubes by integrating over the spectral bands of BLASTPol, 
as follows: 

h S^^(u)T cp (e,u) dv 
Tcp {6) /~E<*M dv • (2) 

Here the superscript "ch" refers to one among 250, 350, and 
500 fim, E ch (v) is the measured spectral response of each of 
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Figure 10. Noise estimation for the spectra shown in Figs.[6]and 
[7] We plot the mean background spectra (in arbitrary units) for 
the co-pol (black solid line) and cross-pol (yellow solid line, shifted 
by 1 in the positive y direction for visual clarity) as a function 
of wavenumber. The (10 it) error bars (in red) are quantified as 
the statistical error on the mean. Also shown for reference is the 
relative spectral response of the three BLASTPol channels, in 
arbitrary units. Henceforth, we adopt a colour code in the plots 
whereby curves referring to the three BLASTPol bands, 250, 350, 
and 500 ^m are drawn in blue, green, and red, respectively. 




50 100 

Rotation angle [deg] 

Figure 11. Band-integrated co-pol modulation functions of the 
BLASTPol HWP at ~120K. The curves show the HWP polari- 
sation modulation functions for a fully polarised source (with a 
fiat spectrum) parallel to the analyser in the three spectral bands. 
Note how the positions of the maxima (and minima) depend on 
the wavelength, even when considering a flat-spectrum polarised 
input source; the dotted vertical lines show the band-integrated 
positions of the HWP extrema (shown in Fig. 1191 1, which result 
from the fitting routine described in the next sections. 



the BLASTPol bands (see lPascale et al.ll2008l ). and T cp (6, v) 
are points on the co-pol surface depicted in Fig. [8] A similar 
expression can be written for the cross-pol band-integrated 
transmission. By performing this integration at every angle 
for which spectral data have been obtained, the interpolation 
of these data points will result in the modulation functions 
of the HWP at ~120K for each of the BLASTPol spectral 
bands; these curves are shown in Fig. [11] for co-pol and in 
Fig- EH f° r cross-pol. 

The modulation curves presented here assume that the 
incoming polarised radiation has no dependence on fre- 
quency, or in other words that the input source has a flat 
spectrum in frequency. Equation [2] can be generalised to in- 
clude the known (or assumed) spectral signature of a given 




40 60 80 100 
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Figure 12. Band-integrated modulation functions equivalent to 
those shown in Fig. Illl but for cross-pol transmission. 



astro nomical or calibration source (see also e.g.. lNovak et al] 
1 19891 . Equation 2). More generally, all the band-averaged 
quantities that we have defined here and will be defined in 
the following are potentially affected by the spectral shape of 
the input source. However, we will see how the HWP trans- 
mission and modulation efficiency are very weakly depen- 
dent on the spectral index of the input source, whereas the 
position of the equivalent axes of the sapphire plate stack is 
m ore significantly aff ected (see also the analysis carried out 
bv lSavini et al.|[2009l ') . especially at 250 and 500 /im. 

Figs. [TT] and [T^] clearly show that there is a significant 
dependence of the position of the HWP maxima and min- 
ima upon frequency, even when considering a flat-spectrum 
polarised input source. These effects are particularly impor- 
tant for a "HWP step and integrate" experiment such as 
BLASTPol, and a careful post-flight polarisation calibration 
must be performed by using all the information available 
from the pre-flight characterisation of the HWP. We begin 
to tackle this problem in the next section, where we outline 
a relatively simple solution to account for most of the HWP 
non-idealities in the data-analysis pipeline, and in particular 
in the map-making algorithm (see Section [5]). 

The spectral transmission datasets of the HWP cooled 
to ~120K, when compared to those taken with the HWP 
at room temperature (Fig. Q, show a definite abatement 
of the in-band losses due to absorption from sapphire, as 
expected. However the effect is still appreciable, especially 
above ~25cm _1 . As we will show in the following, we have 
independent evidence that the residual absorption nearly 
vanishes when the sapphire is further cooled to 4 K, as it is 
when the HWP is installed in the BLASTPol cryostat. While 
it is not currently feasible for us to measure the spectral 
response of the HWP cooled to 4 K, the unique quality and 
completeness of our dataset allow us to fully characterise the 
performance of the BLASTPol HWP. 

We extrapolate our "cold" dataset to 4K, using a com- 
bination of the analytical expression and the data points 
taken from the literature and reported in Fig. [13] The HWP 
modulation efficiency is defined as (T c ° p ° -T°°)/(T° p ° +T°°), 
where the "co-pol" and "cross-pol" transmissions, T c ° p and 
T^p , are the spectral transmission response of the HWP, 
with its axis at 0°, between parallel and perpendicular po- 
larisers, respectively. The inferred co-pol/cross-pol trans- 
missions and modulation efficiency of the BLASTPol HWP 



10 Moncelsi, et al. 



7 0.20 



a 0.15 - 



V, 0.10- 



0.05 




20 25 30 35 40 
Wovenumber [cm -1 ] 



45 50 



Figure 13. Ordinary (solid lines) and extraordinary (dashed 
lines) sapphire absorption coefficient as a function of wavenumber, 
at cryogenic temperatures. The two analytical relations covering 
the whole frequency ra nge are derived by Savini (2010, private 
communication; see also ISavini et aLll2006t) from a set of spectral 
measurements of a sapphire sample at 80 K, and, strictly speak- 
ing, only apply for frequencies < 1 THz (do tted vertical line). 
For re ference, we also plot me asurements fr om Locwcnstci n et all 
l|l973l . diamonds) at 1.5 K and lCook fc Perkowitj l|l985l . squares') 
at 60 K, displaced in x by 0.25 cm" 1 and in y by 0.003 cm" 1 for 
visual clarity; the lines connecting these data points follow the 
convention shown in the legend. 



(with its axis at 0°) at 4K are shown in Fig. 1141 For a 
flat-spectrum input source, the band-integrated transmis- 
sion of the HWP at its maxima is ~0.87, ~0.91, and ~0.95 
at 250, 350, and 500 /im, respectively; whereas the band- 
integrated cross-pol is <0.5%, <0.2%, and <0.5%, respec- 
tively; finally, the band-integrated modulation efficiency is 
-98.8% ~99.5%, and ~99.0%, respectively. 




Wovenumber [cm ] 



(a) The predicted transmissions through the cold HWP as 
a function of wavenumber. The black line shows the HWP 
transmission, T® p , between two parallel polarisers (Q = 1 — ► 
Q = 1) with the HWP axis at 0°. The dark blue line shows 
Q = — 1— > Q = ~ lin the same reference frame (or cquiv- 
alently Q = 1 -> Q = 1 with the HWP axis at 90°). The 
purple line shows the transmission, T® p , with the HWP axis 
at 0° between two perpendicular polarisers. 




Wovenumber [cm '1 
(b) Predicted modulation efficiency of the cold HWP as a 
function of wavenumber, obtained as (T® p — T® p ) / (T® + 
T® p ). Note that the j/-axis scale ranges from 0.8 to 1. 



4 EMPIRICAL MODELLING 

The final goal of this section is to provide a set of usable 
parameters that completely describe the performance of the 
HWP as measured in the laboratory. This set of parameters 
consists of the 16 coefficients of the Mueller matrix of a 
generic HWP, and the actual phase s hift. For an id eal HWP, 
the Mueller matrix at 9 = 0° reads (|CollettJll993T ) 



Mi 



(1 \ 

10 

0-10 

\0 -1/ 



(3) 



and the phase shift is Aip = 180°. 

For a real HWP, these parameters always depart from 
the ideal case to some extent, and certainly depend upon 
frequency. In the following we describe an empirical model 
that we develop specifically for the characterisation of the 
BLASTPol HWP, though we note that it can be applied to 
any HWP to recover its frequency-dependent descriptive pa- 
rameters. Such an empirical model is co mplementary to the 
phys ical and analytical one developed by I Savini et all |2006l . 
120091 ). which produces an analogous output by modelling the 
non- idealities of the components of the HWP assembly and 
their optical parameters. 



Figure 14. Predicted performance of the BLASTPol HWP at 
4K, extrapolated from a set of spectral data collected with the 
HWP cooled to ~120K (see Section 13.311 . "Co-pol" and "cross- 
pol" transmissions, T cp and T xp , are defined as in Fig. [2] The 
approximate extent of the BLASTPol bands is also indicated. 



4.1 Mueller matrix characterisation 

By recalling the Stokes formalism, we can formalise the ex- 
perimental apparatus described in Sections 13.21 and 13.31 as a 
series of matrix products, as follows: 



J out 

QXp 
°OUt 



= D 1 
= D 1 



Ml 
M v p 



TZ(-t 
K(-t 



Mhwp ■ n (9) 
Mhwp ■ Tl (9) 



(4) 
(5) 



Here D is the Stokes vector for a bolometric (polarisation 
insensitive) intensity detector, M p is the Mueller matrix of 
an ideal horizontal polariser, M^ is that of an ideal vertical 
polariser, TZ (9) is the generic Mueller rotation matrix, and 
Sin is the horizontally polarised input beam from the pFTS. 
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By expanding all the matrices in Equation [4] 
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we can compute the products and rearrange as 



Sout = 2 [ 2a °° + ai1 + a22 + 2 ( a<)1 + ai °) cos 2 " + 
+(an — 022) cos 46" — 2(ao2 + 020) sin 28 — (ai2 + 021) sin4#j 
= A + Bsin26» + Ccos26l + Dsin46l + £;cos46l , 

with 

A = a 00 + — + — , 

B = — (a 2 + a 2 o) , C = aoi + aio , (9) 
D = — -(012+021), S=i(an — 022) ■ 
Similarly, noting that 



-Mp = 



/ 1 


-1 





°\ 


-1 


1 




























0/ 



we rearrange Equation [5] as 

Soft = A' + B' sin 26» + C" cos 20 + L>' sin 46* + E' cos 46* , 
with 

./ Ill 122 

= aoo- — -—, 

B = a2o — 002, C = aoi — «io , (11) 

D' = -(012 + 021), E' = — (022 — an) . 

Finally, by performing linear combinations of the quan- 
tities defined in Equations [5] and 1111 one can write the 
individual elements that compose the Mueller matrix of a 
generic HWP as 

aoo = \ {A + A') , aoi = | (C + C') , 
aio = i(C-C'), 011 = i (A-A'+E-E') , (12) 
a 02 = -~ (B + B'), a 20 = i {B 1 - B) , 

a 22 = ~ (A- A' - E + E') , oia = a 2 i = \ ip' - D) , 

where in the last equality we currently assume the symmetry 
of two coefficients, 012 = 021. This degeneracy may be bro- 
ken by imposing the conservation of energy, i.e. by requiring 
the output Stokes vector resulting from a generic polarised 
input travelling through the recovered HWP Mueller ma- 
trix to satisfy I 2 = Q 2 + U 2 . Alternatively, the degeneracy 



can be broken by taking spectra at an intermediate configu- 
ration between co- and cross-pol; this additional constraint 
will be included in a future work (Spencer et al., in prepa- 
(6) ration). Also, because our experimental setup is sensitive to 
linear but not circular polarisation, this method only allows 
us to constrain the 9 elements of the Mueller matrix associ- 
ated with [I,Q,U]. The remaining 7 coefficients associated 
with V can only be measured with the use of a quarter- 
wave plate, which induces a phase shift of 90° between the 
two orthogonal polarisations travelling through the plate; 
this measurement is beyond the scope of this paper and not 
pertinent to the needs of BLASTPol. 

We want to estimate the 9 coefficients derived in Equa- 
tion [12] from the co-pol and cross-pol data cubes described 
in Section 13.31 Equations [5] and [TU] encode a simple depen- 
dence of an d upon 6, the HWP rotation angle. 
Therefore, for a given frequency, a fitting routine can be 
' applied to the measured transmission curves as a function 
of 8, to determine the parameter sets [A, B, C, D, E] and 
[A', B' , C' , D' , E'] for the co-pol and cross-pol configura- 
(8) tions, respectively. By repeating the fit for every frequency, 
we have an estimate of the 9 coefficients as a function of 
wavelength. However, this procedure does not allow us to 
associate any uncertainty to our estimates. 

A better approach to this problem is to use a Monte 
Carlo simulation. We repeat the above fitting procedure 
1000 times; every time we add to every individual trans- 
mission curve a realisation of white noise, scaled to the 1 a 
spectral uncertainty as estimated in Fig. 1101 and compute 
the fit using this newly generated transmission curve. In ad- 
dition, for every frequency we introduce a random jitter on 
the rotation angle that has a la amplitude of 1°. The dis- 
persion in the fitted parameters due to these two types of 
uncertainties, which are inherent to the measurement pro- 
cess, provides a realistic estimate of the uncertainty asso- 
ciated with each of the 9 coefficients. In particular, at each 
'frequency, we produce 9 histograms of the 1000 fitted values. 
We use the mode of each distribution as our best estimate 
for the corresponding coefficient at that frequency, and the 
68% confidence interval as the associated 1 a error. 

In Fig. [TS] we show a graphical representation of the 9- 
element Mueller matrix of the BLASTPol HWP at a given 
angle (6* = 0°), as a function of wavenumber. In Figs. 1 1 611 1 Tt 
and[TH] we show the resulting histograms for the 9 coefficients 
at 20, 28.6, and 40 cm -1 , respectively (which are the central 
frequencies of the BLASTPol bands). 

4.2 Position of the HWP equivalent axes 

The behaviour of the coefficients as a function of wavenum- 
ber shown in Fig. [T5] suggests that the position of the HWP 
equivalent axes, /3 C a hereafter, may have an inherent fre- 
quency dependence, which we must investigate. /3 ea can be 
readily retrieved at each frequency by locating the rotation 
angle that corresponds to the first minimum in the fitted 
transmission curve. Hence, /3 oa is measured with respect to 
an arbitrary constant offset that is inherent to the specific 
experimental setup; we set this offset to be zero at 25 cm -1 . 
Operationally, this means that the HWP zero angle in the 
instrument reference frame (/3o; see Equation I17p must be 
calibrated using the 350 /im band. A plot of j3 ea , as a function 
of wavenumber is given in Fig. 1191 
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Figure 15. Graphical representation of the Mueller matrix of 
the BLASTPol HWP at a given angle (6 = 0°), as a function of 
wavenumber. The (10 cr) error bars (in red) are quantified via a 
Monte Carlo, which accounts for random errors in the spectra of 
amplitude as given in Fig, 1101 and random errors of amplitude 1° 
in the rotation angle. 



As anticipated, it is of crucial importance to derive the 
band-averaged value of /3 ea for input sources with different 
spectral signature, as follows: 



P ™ J °° E ch (u) q (y) dv 



(13) 



where we adopt the same notation as in Equation [2] and 
the known (or assumed) spectrum of an astronomical or 
calibration source is modelled as c (y) oc v a . We compute 
Equation[l3]for a range of spectral indices of interest: a = 
for a flat spectrum; a = 2 for the Rayleigh- Jeans tail of a 
blackbody; a — 4 for interstellar dus t, modelled as a m odi- 
fied blackbody with emissivity /3 = 2 ijHildebrandl 19831 ) ; and 
finally a — — 2 as a replacement for the mid-infrared expo- 
nential on the Wien side of a blackbody to acc ount for the 
variability of dust temperatures within a galaxy dBlainll 19991 : 
iBlain et al.ll2003t ). The results of this analysis are shown in 
Fig. QI] and in Tabled] 

As expected, the impact of different input spectral sig- 
natures is minimal at 350 /im, where the HWP has been de- 
signed to function optimally (see Section [2.2[) : whereas the 
spectral dependence is more pronounced at 250 and 500 /im, 
and, if neglected, it may lead to an arbitrary rotation of 
the retrieved polarisation angle on the sky of magnitude 
2/3 ca = 10-15° (3-5°) at 250 (500) /im (see Equation . 

We have thus confirmed that the dependence of the 
HWP equivalent axes upon wavelength is inherent to the 
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Figure 16. Histograms at 20 cm 1 (central frequency of the 
500 /xm BLASTPol band) resulting from the Monte Carlo fit of 
the HWP parameters. For every histogram, the dashed red line 
indicates the mode of the distribution, which we adopt as our 
best estimate for the corresponding coefficient at that frequency, 
while the two dotted red lines indicate the 68% confidence inter- 
val, which we use as the uncertainty on the retrieved coefficient. 
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Figure 17. Equivalent histograms to those shown in Fig. |16l but 
at 28.57cm -1 (central frequency of the 350 fim BLASTPol band). 
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Figure 18. Equivalent histograms to those shown in Fig. I16l but 
at 40 cm - 1 (central frequency of the 350 /an BLASTPol band). 
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Figure 19. Position of the HWP equivalent axis, /3 ea , as a func- 
tion of wavenumber (solid black line) . Note that this quantity is 
defined with respect to an arbitrary constant offset that is inher- 
ent to the specific experimental setup; we set this offset to be zero 
at 25 cm - 1 . The band-averaged values for input sources with dif- 
ferent spectral index (a; see legend) are drawn as thick horizontal 
lines. Also shown for reference is the relative spectral response of 
the three BLASTPol channels, in arbitrary units. 
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0.44 


1.6 



Table 1. Band-averaged position of the HWP equivalent axis for 
sources with different spectral index. The input source is assumed 
to have a spectrum ? oc v a . 
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Figure 20. Graphical representation of the Mueller matrix of 
the BLASTPol HWP, equivalent to that shown in Fig. [T5j but 
including in the fit the frequency-dependent position of the HWP 
equivalent axes (see Fig. 1 19D . 



achromatic design. We now postulate that most of the non- 
idealities we see in the measured HWP Mueller matrix 
(Fig. I15[l are primarily due to the wavelength dependence 
of pea., along with the residual absorption from sapphire at 
~120K. This hypothesis naturally ensues from the discus- 
sion presented in Section 12.11 on the scatter in frequency 
that results from any polarisation rotation on the PS sphere 
produced by a multiple-slab wave plate. The measurements 
of pea. presented in Fig. [19] effectively quantify the area of 
the PS surface in which the various polarisation states re- 
group. One can imagine that the HWP performance would 
approach the ideal case once this effect is corrected for. 

Therefore, we include Pea (v) in our Monte Carlo as a 
frequency-dependent offset in the array of rotation angles 
(so that 9-^9 — Pea), and repeat our simulations. The re- 
sults, presented in Fig. 1201 can now be qualitatively com- 
pared to the Mueller matrix of an ideal HWP (Equation [3]). 
The improvement is noticeable, especially in the off-diagonal 
elements, and the resemblance to an ideal HWP is remark- 
able across the entire spectral range of interest; this proce- 
dure effectively acts to diagonalise the HWP Mueller matrix. 
However, the transmission losses due to absorption from the 
sapphire at ~120K still affect the diagonal elements of the 
matrix, as expected. 

As a final improvement, we extrapolate the p ca - 
corrected HWP Mueller matrix to 4 K by including in our 
Monte Carlo a correction for the residual sapphire absorp- 
tion (using the data presented in Fig. I13[) . The results are 
shown in Fig. [21] Although there still seems to be residual 
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Figure 21. Graphical representation of the Mueller matrix of the 
BLASTPol HWP, equivalent to that shown in Fig. [20] but includ- 
ing a correction for the temperature dependence of the sapphire 
absorption coefficient (see Fig. 1131 1. 



transmission losses due to sapphire absorption at 250 and 
350 fim, the retrieved HWP Mueller matrix is nearly that 
of an ideal HWP. The band-averaged values of the matrix 
coefficients for a fiat-spectrum input source are reported in 
Table [21 along with their propagated uncertainty; the off- 
diagonal elements are always consistent with zero within 2 a 
and the modulus of the three diagonal coefficients is always 
> 0.8. The combination of these coefficients with the band- 
averaged values of /3 ca given in Table [1] gives a complete 
account of the HWP non-idealities to the best of our ability. 

We repeat the calculation of the band-averaged coeffi- 
cients for the other spectral indices discussed in Fig. 1191 we 
find values that are always within 1-2% of those reported in 
Table [21 and thus we do not explicitly report them here. Be- 
cause the three diagonal elements of the HWP Mueller ma- 
trix effectively determine the HWP co-pol/cross-pol trans- 
mission and modulation efficiency, this analysis confirms 
that these quantities are very weakly dependent on the spec- 
tral index of the input so urce; these findings are in very good 
agreement with those of ISavini et al l (|2009h . We will SCO 111 
Section [5] how Soo, flu, and 7I22 can be incorporated in the 
map-making algorithm in terms of optical efficiency, rj, and 
polarisation efficiency, e, of each detector. 



4.3 Effective HWP phase shift 

Finally, we discuss a potential limitation to any linear polar- 
isation modulator, i.e. the leakage between axes. In a HWP, 
the phase shift between the two axes should be as close to 



Band 


250 £im 


350 /Jin 


500 [im 


«oo 


0.905 ± 0.006 


1.001 ± 0.006 


1.008 ± 0.007 


aoi 


0.012 ± 0.010 


0.017 ± 0.010 


0.014 ± 0.011 


Q02 


-0.002 ± 0.008 


0.006 ± 0.009 


0.001 ± 0.009 


aio 


-0.016 ± 0.010 


-0.021 ± 0.010 


-0.020 ± 0.011 


an 


0.806 ± 0.011 


0.928 ± 0.010 


0.935 ± 0.012 


a 12 


-0.007 ± 0.011 


-0.009 ± 0.014 


-0.011 ± 0.014 


Q20 


-0.008 ± 0.008 


-0.022 ± 0.010 


-0.021 ± 0.010 


021 


-0.007 ± 0.011 


-0.009 ± 0.014 


-0.011 ± 0.014 


Q22 


-0.808 ± 0.008 


-0.960 ± 0.009 


-0.979 ± 0.010 



Table 2. Band-averaged Mueller matrix coefficients. These values 
are relative to Fig. 1211 The input source is assumed to have a flat 
spectrum in frequency. 



180° as possible to avoid transforming linear polarisation 
into elliptical polarisation, hence losing efficiency. The phase 
can not be directly measured in a pFTS, but it can be indi- 
rectly inferred from the HWP Mueller matrix. 

In order to recover the wavelength-dependent phase 
shift of the HWP, we recall the Mueller matrix of a non-ideal 
impe dance-matched single birefringent slab (jSavini et alJ 
2009, at 9 = 0°): 



0' 



',A^) = -x 



(14) 



(a 2 + p 2 
a 2 -? 2 


V 



a 



2+/3 2 









2 a j3 cos Aip 
—2 a /3 sin Acp 



2 a j3 sin Aip 
2a P cos Aip ) 



By comparing the matrix in Equation [14] with that of a 
generic HWP, we can solve for the HWP phase shift as 



cos Atp = 

r 2 



^ app + api j 2 ^ apo — oqi j ' 



2 J V 2 J ™ 
Equation [15] allows us to recover the phase shift from 
our knowledge of app, ap\ and 022- Fig. 1221 shows the esti- 
mated phase shift of the BLASTPol HWP as a function of 
wavenumber, before and after the introduction in our Monte 
Carlo routine of the wavelength-dependent position of the 
HWP equivalent axes depicted in Fig. 1191 The improvement 
is striking, and confirms the fact that most of the HWP non- 
idealities due to the achromatic design can be more easily 
modelled by estimating /3 ca (v) . This finding further encour- 
ages us to implement /3 ca in the map-making code. 

Nonetheless, the /3 ea -corrected phase shift appreciably 
departs from 180°. We have anticipated that this discrep- 
ancy is primarily due to the ~0.047mm difference between 
the desired thickness of the single sapphire substrates and 
what was available on the market; at 350 ^m, this trans- 
lates to a deviation of ~15° from ideality, which accounts 
for roughly 60-75% of the measured non-ideality. The re- 
maining 25-40% can be qualitatively explained by recalling 
that Equation[l5]strictly applies only to a single birefringent 
plate. In a multi-slab Pancharatnam stack, Aip becomes an 
"effective" phase shift as we are no longer in the presence 
of a pure cosine modulation, thus slightly skewing the HWP 
modulation function and resulting in an artificially higher 
leakage between axes when the cosine function is inverted. 

However, we have indications that the modulation effi- 
ciency of the HWP at 4 K is only mildly affected by this 
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Figure 22. HWP phase shift as a function of wavenumber, be- 
fore (orange) and after (black) implementing in the Monte Carlo 
the wavelength-dependent position of the HWP equivalent axes 
(Fig, 119b . The (3<r) error bars (in yellow) are obtained by propa- 
gating the error on the Mueller matrix coefficients. The band- 
averaged values of the phase shift (for a flat-spectrum input 
source) are drawn as thick horizontal lines (only for the upper 
black line). Also shown for reference is the relative spectral re- 
sponse of the three BLASTPol channels, in arbitrary units. 

departure from ideality. From Fig. 114b we see that the 
extrapolated HWP modulation efficiency is always above 
95% across the whole spectral range of interest, with band- 
integrated values exceeding 98%. Moreover, phase shift de- 
viations of similar amplitude are measured in most mm and 
submm-w ave achromatic ha l f- wave plates man ufactured to 
date (e.g.. lSavini et al.ll2009l : IZhang et al.ll201ll ) 

Incidentally, we verify that our methodology does not 
violate conservation of energy by ensuring that the output 
Stokes vector resulting from a generic polarised input trav- 
elling through the recovered HWP Mueller matrix satisfies 
I 2 ^ Q 2 + U 2 in every instance describe above. 



5 MAP-MAKING ALGORITHM 

Map-making is the operation that generates an astronom- 
ical map, which contains in every pixel an estimate of the 
sky emission, and is obtained by combining data from all de- 
tectors available at a given wavelength channel, their noise 
properties and the pointing information. The raw data con- 
sist of bolometer time-ordered streams (or timelines), which 
are cleaned and pre-processed before being fed into the map- 
maker: in order, cosmic rays are flagged and removed, the 
known electronics transfer function is deconvolved from the 
data streams, an elevation-dependent common-mode sig- 
nal due to the residual atmosphere is removed concurrently 
with a polynomial fit to the data, and finally the timelines 
are high-pass filtered to suppress the low-frequency (1//) 
noise. The details of the pre-processing of the BLAST time- 
lines are extensively described elsewh ere (|RexH2007l ; iTruchl 
l2007l ; IWiebell2008l ; |Pascale et al.|[20oj ). and we refer to these 
works for a complete account of the low-level data reduction. 
The process of cleaning and preparing the bolometer time- 
streams for map-making in BLASTPol has closely followed 
that of BLAST, except for the removal of discontinuities in 
the DC level of the bolometer, caused by the HWP being 
stepped approximately every 15 minutes (this operation is 



performed before the high-pass filtering); also, the subtrac- 
tion of an elevation-dependent term from the timelines was 
not needed in BLAST. 

In the following, we focus on the mathematical formal- 
ism of the map-making technique, and its algorithmic im- 
plementation in the specific case of BLASTPol. 

5.1 Maximum likelihood map-making 

For a non-ideal polarisation experiment, by adopting the 
Stokes formalism and assuming that no circular (V) polari- 
sation is present, we can model the data as 

4 = ^Al p [l p +e l (Q p cos2 7t ' + U p sin2 7 ;)] + n\ . (16) 

Here: i, t and p label detector index, time, and map pixel re- 
spectively; d\ are the time-ordered data for a given channel, 
related to the sky maps [I p , Q p , U p ] by the pointing operator 
A\ p \ rf is the optical efficiency of each detector; e l is the po- 
larisation efficiency of each detector with its polarising grid 
(analyser); and n\ represents a generic time-dependent noise 
term. Throughout this discussion it is assumed that the term 
within square brackets is the convolution of the sky emis- 
sion with the telescope point-spread function (PSF). 7^ is 
the time-ordered vector of the observed polarisation angle, 
defined as the angle between the polarisation reference vec- 
tor at the sky pixel p (in the chosen celestial frame) and the 
polarimeter transmission axis. 7I is given by 

ll = al + 2 [p t -p -p J +5' gTid , (17) 

where a\ is the angle between the reference vector at pixel 
p and a vector pointing from p to the zenith along a great 
circle, fit is the HWP orientation angle in the instrument 
frame, /3o is the HWP zero angle in the instrument frame, 
Pes. ls t ne band-averaged position of the equivalent axes of 
the HWP (dependent on the known or assumed spectral 
signature of the input source; see Section [4.2|) . and <5g rid = 
[0, 7r/2] accounts for the transmission axis of the polarising 
grids being parallel/perpendicular to the zenith angle. 

The notation outlined above can be connected to the 
Mueller formalism developed in Section R.ll to determine un- 
der which circumstances Equationfjniis valid in the presence 
of a real (i.e., non-ideal) HWP. Because we have included in 
Equation [17] the band-averaged position of the equivalent 
axes of the HWP, /3 ca , the Mueller matrix of the BLASTPol 
HWP can be considered almost that of an ideal HWP, as 
discussed in Section 14.21 Nonetheless, we have shown that 
the band-averaged values of the three diagonal matrix coef- 
ficients are not identically unity (but always > 0.8 in mod- 
ulus), probably as a result of residual absorption from sap- 
phire, especially in the 250 and 350 /im bands (although we 
have corrected for it to the best of our knowledge). 

In the light of these considerations, we now want to com- 
pare Equation [TB] to Equation [SI which both represent the 
signal measured by a polarisation insensitive intensity detec- 
tor when illuminated by a polarised input that propagates 
through a rotating HWP and an analyser. A term-by-term 
comparison shows that these two expressions are equivalent 
when the coefficients B and C (defined in Equation [9]) are 
zero, i.e. when the HWP modulates the polarisation purely 
at four times the rotation angle, with no leakage in the sec- 
ond harmonic (twice the rotation angle) and thus no leakage 
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of I into Q and U. These two coefficients are linear combina- 
tions of the HWP Mueller matrix elements aoi, aio, 0,02, 020, 
which we have shown in Table [2] to be all compatible with 
zero within 2 a. In addition, their amplitude is at most ~2% 
of that of the diagonal matrix elements, and in the limit of el- 
evated polarisation angle coverage, (cos2 7 ) 2 + (sin2 7 ) 2 « 0, 
these terms (in twice the rotation angle) effectively average 
out in the sums. Therefore, the coefficients B and C can 
be neglected to first order, and the two expressions can be 
considered equivalent. Nonetheless, these generally moder- 
ate levels of I — > Q, U leakage can be readily accounted for 
by incorporating in the map-making algorithm a correction 
for the "instrumental polarisation" (IP). 



In addition, after some elementary algebra, it can be 
shown that 77 = a 00 + ^ + ^f - and ne = 211 _ 222. As an- 
ticipated in Section r4.2l the knowledge of the band-averaged 
values of the three diagonal matrix elements, aoo,«i 1,022 
(which we have shown to depend weakly on the spectral in- 
dex of the input source), can be readily incorporated in the 
map-making algorithm in terms of optical efficiency, n, and 
polarisation efficiency, e, of the HWP; these can be factored 
into the overall optical efficiency and polarisation efficiency 
of each detector. From the values listed in Table [2] in our 
case we find [??hwp,£hwp] = [0.904,0.893], [0.985,0.958], and 
[0.986, 0.971] at 250, 350, and 500 fim, respectively. 



Finally, the comparison of Equations [16] and Equation[8] 
also yields nex = —112 = —021, where we have introduced 
a new parameter, x> which quantifies the amplitude of the 
mixing of Q and U. From Table O we see that 012 = «2i 
are always compatible with zero within I a, and their am- 
plitude is at most ~1% of that of the diagonal matrix ele- 
ments. Nonetheless we quantify the amplitude of the Q <H> U 
mixing to be Xhw P = 0.009,0.010, and 0.011 at 250, 350, 
and 500 /im, respectively. While this correction is not cur- 
rently included in our algorithm, we indicate that it can be 
implemented in a relatively straightforward way by mod- 
ifying Equation [16] with a double change of variable, i.e. 
Q — s> Q + \U and U — > U + xQ- If X ls estimated to the 
required accuracy, the unmixed Q and U can be retrieved 
unbiasedly. This correction may be very relevant to CMB 
polarisation experiments, where any Q o U leakage leads 
to a spurious mixing of the EE and BB modes. 



We remind the reader that the above factors have been 
computed directly from the band-averaged coefficients of the 
inferred HWP Mueller matrix extrapolated to 4K, and of- 
fer a direct way to include the modelled HWP non-idealities 
in a map-making algorithm. On the other hand, the band- 
averaged HWP maximum transmission, polarisation effi- 
ciency and cross-pol quoted at the end of Section 13.31 are 
estimated directly from the spectra extrapolated to 4K, and 
are mostly informative from an experimental point of view 
rather than for data analysis purposes. 



Consider now one map pixel p that is observed in one 
band by k detectors [i — l,...,fc); let us define the gener- 
alised pointing matrix A tp , which includes the trigonometric 



functions along with the efficiencies, 

/V Al p n 1 e 1 Alp cos 2-yl n 1 e 1 A\ p sin 2 7t a \ 



Hp 



Hp 



rf e { A\ v cos 2-yl rf e l A\ v sin 2 7 



\n h A% V k e k A* p cos2 7t fc r, k s k A% sin 2 7t fc / 

(18) 

and the map triplet S p , along with the combined detector 
(T>t) and noise (nt) timelines, 



fdl\ 
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\d h ) 



(nl\ 



(19) 



Equation [16] can then be rewritten in compact form as 

V t = A tp S p + n t . (20) 

Under the assumption that the noise is Gaussian and 
stationary, the likelihood of S p given the data can be 
maximised, thus yielding the well-known generalised least 
squares (GLS) estimator for S p : 



AZpN^A,,,) Aj v 1SF 1 Vt 



(21) 



where N is the noise covariance matrix of the data in the 
time domain, 



N= (n t rv> = 



((nUl,) 



{nln],) 



\{n k t nl, 



(nl nl, 



{n\n\,) 



(n k nl,} 



(nln k ,)\ 



(nln k ) 



(n k n k ,)J 



(22) 

with t, t' running over the detector time samples (typically 
iV s ~ 10 B -10 7 ). 

Computation of the solution to Equation [21] is far from 
trivial in most astronomical applications, due to N being a 
very large matrix, of size kN s x kN s . Understandably, it is 
computationally challenging to invert this matrix, especially 
when there are correlations among detectors, and a number 
of "optimal" map-making techniques have be en developed 



in the literature to tackle this problem (e.g., Natoli et al 
l200ll. l2009l:lMasi et al.ll2006l: Ijohnson et"aLjl2007l: IWu et al 



120071 : IPatanchon et a l. 2008; Canta lupo et al.ll2010h 



5.2 Naive binning 

In the simple case that the noise is uncorrelated between 
different detectors, the matrix in Equation (|22[) reduces to 
block diagonal: 



{n\n 3 t ,) = (n J t n\,) = (i / j) 



(23) 



In addition, let us assume that there is no correlation be- 
tween noise of different samples acquired by the same de- 
tector, or, in other words, that the noise in each detector 
is white. From Equations (|22|) and (I23[l . we can see that 
each "block" of the noise covariance matrix collapses into 
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one value, which is the timeline variance for each detector. 
Hence, N becomes a k x k diagonal matrix where the di- 
agonal elements are the sample variances of the detectors, 
of, and weights can thus be defined as the inverse of those 
variances, w t = 1/of. 

Therefore, under the assumption that the noise is white 
and uncorrelated among detectors, Equation (|21[1 reduces 
to a simple, weigh ted binning ("naive" binning; see also 



iPascale et aT]|201ll ) of the map: 




E w i i 
i=it=i 1 tp' *p 

k 

E wi 
i=i 



(24) 



If we now define the quantities, 

N hi t = 2' c = E \ cos 27 *' 

i,t i,t 

C2 = \ cos2 2 7t, s = Yl \ sin 27 *' 

i,t i,t 

S2 s ^2 o sin2 2 T* = ^ hit ~ C2 ' U = ~Y, u t, 

i,t i,t 

m = ^ - cos 27] sin 27! , = ^ Ut cos 27? 

S2 = Y Ut sm2 Tt> 



C2 =Y1 fP* cos 2 7* ' S 2 = £ ^ sin 2 T* ' 



(27) 



In the light of these considerations, let us go back to 
Equation (|16|l and model the generic time-dependent noise 
term nl as 



n' t =u t + Cpt , 



(25) 



where ut represents a time-dependent noise term, completely 
uncorrelated among different detectors, while pt describes 
the correlated noise (varying over timescales larger than the 
ratio of the size of the detector array to the scan speed), 
coupled to each detector via the £ l parameter, peculiar to 
each bolometer. 

Let us define the following quantity for every pixel p in 
the map: 




1 E E dl \ 

i=l t = l 

k N s 

E E rft cos2 7 j 

! = 1 t = l 

k N s 



y E E sin2 7 ; 

\i=l i=l 



(26) 



where N s is now the number of samples in each detector 
timeline that fall within pixel p, and the superscript "e" 
stands for "estimated". The above quantities can be com- 
puted directly from the detector timelines. Recalling Equa- 
tions (|I6p and (|25p . we can outline the following linear sys- 
tem of 3 equations with 3 unknowns: 



1 US j 
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where we have temporarily assumed rf = e x = w 1 — 1 and 
combined the two sums in one, with the indices i and t 
running, respectively, over the bolometers and the samples 
in each detector timeline. 



then the system in Equation (|27l) can be rewritten in com- 
pact form as 
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(28) 

In order to retrieve an estimate of S p from the quan- 
tities computed in Equation (I26[) . the above system has to 
be solved for every pixel p in the map. One can already 
see the computational advantage of inverting a 3 x 3 matrix 
A^pix x TVpix times, with respect the inversion of a generic 
kN s x kN a matrix (for detectors having uncorrelated l/,f 
noise as well as a common-mode 1// noise; IPatanchon et all 
2008), or k matrices of size N s x N a (for detectors having 
only uncorrelated 1// noise: ICantalupo et al.l [2010f). The 
main difficulties are, of course, in estimating the noise terms 
U,P,C%,C%,S%,S%. However, recalling Equation {T7J) and 
the fact that adjacent detectors have orthogonal polarising 
grids (<5g r i d = [0, t/2]), we note that, in the sum over i, adja- 
cent detectors have equal and opposite contributions to 
and Sji under the following assumptions: (i) the timescale 
over which the correlated noise is approximately constant is 
larger than the time elapsed while scanning the same patch 
of sky with two adjacent detectors, and (ii) f; 1 is not too 
dissimilar between adjacent bolometers. 



This means that the terms Cj and 



5*2 can be neglected, 



under the above assumptions, while estimating the [Q, U] 
maps. In particular, as a first step, we can solve for / only by 
high-pass filtering the timelines, in order to suppress the cor- 
related noise term in /, P. Subsequently, / can be assumed 
known, and the [Q, U] maps can be computed without filter- 
ing the timelines, so that polarised signal at large angular 
scales is not suppressed. In fact, we see from Equation (|16[1 
that in the limit of elevated angle coverage, the term in I, 
not being modulated at four times the HWP rotation angle, 
effectively averages out in the sums. 

The other assumption required for the naive binning is 
that the noise is white, at least on the timescales relevant 
to BLASTPol's scan strategy. An analysis of the bolometer 
timelines from the 20f campaign shows that the knee of the 
1/f noise in the difference between two adjacent detectors is 
typically located at frequencies < 0.1 Hz; assuming a typical 
scan speed of 0.1° s _1 , this corresponds to angular scales of 
> 1° on the sky. The regions ma pped by BLAS TPol hardly 
exceed I ° in size (see Table 1.1 in 1]VIoncelsill201lf ). hence here 
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we stipulate that the noise in the difference between pairs 
of adjacent detectors is white. 

Therefore, under the assumptions above, we can solve 
the linear system outlined in Equation (|28f) : by defining the 
following quantities, 

A = c (c 2 - N h ) - N h (cl + m 2 — c 2 N h ) + 2csm ~ c 2 s 2 , 

A = - (c 2 + m 2 - c 2 Nh) , B = c (c 2 — Nh) + sm , 

C = cm~sc 2 , D = — [(c 2 — Nh) Nh + s 2 ] , (29) 

E = cs-mNh, F = c 2 Nh~c, 

the solution to the system can be written in compact form: 

/ r X / AIl+BQ° p +CU? 



Q P = 



BI° P +DQ° P +EU^ 



(30) 



where we have renamed iVhit —> Nh for brevity. 
5.3 Weights and uncertainties 

The solution for S p given in Equation (|30[) is a simple, un- 
weighted binning of the data into the map pixels. In real- 
ity, as anticipated in Equation (|24[) . we want to perform a 
weighted binning, where the weight of each detector is given 
by the inverse of its timeline variance, which can be easily 
measured as the bolometer's white noise floor level. In our 
formalism, the weighted binning is simply achieved by defin- 
ing [Ip, Q P , Up] in Equation (|26p . as well as each of the quan- 
tities Nh, c, s, c 2 , S2, and m introduced in Equation (|27[) . 
to include w l in the sums. Similarly, the measured values 
of the optical efficiencies rf and polarisation efficiencies e* 
can readily be inserted in Equation (|27[l to account for the 
non-idealities of the optical system. 

The introduction of the weights allows us to derive the 
expression for the statistical error on S p , in the continued 
assumption of uncorrelate d noise, following the usual error 
propagation formula (e.g., iPress et al, I ll992l . here we omit 
the sum over t for simplicity): 



(31) 



After some tedious algebra, the expression for the statistical 
error is 

/Vai£\ 
Op = Varp 3 = 
\Varp7 

/■A (A 2 N h + B 2 c 2 + C 2 s 2 + 2ABc + 2AC s + 2BCm) 
(B 2 N h + D 2 c 2 + E 2 s 2 + 2 B D c + 2 B E s + 2 D E m) 
V % (C 2 Nh + E 2 c 2 + F 2 s 2 + 2 C E c + 2 C F s + 2 E F m) 

where S2 = Nh—Q2, as noted in Equation (|27|) . To first order, 
these expressions can be used to quantify the uncertainty of 
[I, Q,U] in each map pixel p. A more comprehensive account 
of the correlations in the noise, as well as a thorough valida- 
tion of the assumptions made here, is beyond the scope of 
this paper and will be treated in a future work. 



5.4 Preliminary map 

As a proof of concept of the naive binning technique for the 
BLASTPol polarised map-maker, which includes all the cor- 
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Figure 23. Preliminary BLASTPol map at 500 ^m of the Ca- 
rina Nebula, a giant molecular cloud approximately centred at 
coordinates [10 h 42 m 35 s , -59°42'15"]. The intensity contour lev- 
els are shown in the background. The pseudo- vectors indicate the 
inferred magnetic field direction: the BLASTPol measurements 
(red) are in excellent agree ment with thos e from the SPARO po- 
larimeter at 450 fim (black; iLi et alj [2006V This map should not 
be considered of any scientific value as no absolute calibration has 
been applied to the intensity contours; the map is only shown as 
a proof of concept for the map-maker, which includes all the cor- 
rections due to HWP non-idealities discussed in this work. Here 
we adopt /3 ea ^ m =1.7, cal culated assuming a s pectral index for 
the dust emissivity of 1.37 llSalatino et al.l 2012 V 



rections due to HWP non-idealities discussed in this work, 
we present a preliminary map at 500 fim of one of the giant 
molecular clouds observed by BLASTPol, the Carina Neb- 
ula. In this specific case, we calculate /? ea M m = 1.7 based 
on a dust emissivity spectral index of 1.37 (|Salatino et al.l 
2012) for a modified blackbody spectral energy distribution. 

The map in Fig. [23] is presented as contour levels of 
the intensity map /, upon which we superimpose pseudo- 
vectors indicating the inferred magnetic field direction on 
the sky (assumed to b e perpendicula r to the measured po- 
larisation direction; see lLazarianll2007l 'l . The sky polarisation 
angle is given by <j) — i arctan ^ . Because the absolute flux 
calibration has not been finalised yet, we choose not to re- 
port here the intensity values corresponding to each contour 
level. This map should not be considered of any scientific 
value as it is not calibrated in flux. Nonetheless, we find 
a remarkable agreement between the BLASTPol polarisa- 
tion directions at 500 /im and those produced by the Sub- 
millimete r Polarimeter for Antarctic Remote Observations 
(SPARO; iNovak et al.ll2003h at 450 rim, which are shown in 
Fig. 1 of ILi et al.l l|2006r ) and whose pseudo- vectors are re- 
ported in Fig. [23] for comparison. Here the original BLAST- 
Pol [I,Q,U] maps have been smoothed with a kernel of 4' 
(FWHM) for a more direct comparison with the SPARO 
data. 



6 CONCLUDING REMARKS 

The goal of the first part of this work was to identify and 
measure the parameters that fully characterise the spectral 
performance of the linear polarisation modulator integrated 
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in the BLASTPol instrument, a cryogenic achromatic half- 
wave plate. 

We have described in detail the design and manufac- 
turing process of a five-slab sapphire HWP, which is, to our 
knowledge, the most achromatic built to date at submm 
wavelengths. In the same context, we have provided a useful 
collection of spectral data from the literature for the sap- 
phire absorption coefficient at cryogenic temperatures. 

Using a polarising FTS, we have fully characterised the 
spectral response of the anti-reflection coated BLASTPol 
HWP at room temperature and, unprecedentedly, at 120 K; 
we have acquired data cubes by measuring spectra while 
rotating the HWP to produce the polarisation modulation. 

The cold dataset contains measurements in both co-pol 
and cross-pol configurations; we have used these two data 
cubes to estimate 9 out of 16 elements of the HWP Mueller 
matrix as a function of frequency. We have developed an ad- 
hoc Monte Carlo algorithm that returns for every frequency 
the best estimate of each matrix element and the associ- 
ated error, which is a combination of the uncertainty on the 
measured spectra and a random jitter on the rotation angle. 

We have measured how the position of the equivalent 
axes of the HWP, /3 e a, changes as a function of frequency, an 
effect that is inherent to any achromatic design. Once this 
dependence is accounted for in the Monte Carlo, and a cor- 
rection is implemented for the residual absorption from sap- 
phire, the Mueller matrix of the HWP approaches that of an 
ideal HWP, at all wavelengths of interest. In particular, the 
(band-averaged) off-diagonal elements are always consistent 
with zero within 2 a and the modulus of the three diagonal 
coefficients is always > 0.8. Therefore, we have introduced in 
the BLASTPol map-making algorithm the band-integrated 
values of ,9 ca as an additional parameter in the evaluation of 
the polarisation angle. To first order, this approach allows 
us to account for most of the non-idealities in the HWP. 

We have investigated the impact of input sources with 
different spectral signatures on /3 ca and on the HWP Mueller 
matrix coefficients. We find that the HWP transmission 
and modulation efficiency are very weakly dependent on 
the spectral index of the input source, whereas the position 
of the equivalent axes of the sapphire plate stack is more 
significantly affected. This latter dependence, if neglected, 
may lead to an arbitrary rotation of the retrieved polarisa- 
tion angle on the sky of magnitude 2/3 ca = 10-15° (3-5°) 
at 250 (500) /im. The 350 /im band, however, is minimally 
perturbed by this effect. 

In principle, the measured Mueller matrix can be used 
to generate a synthetic time-ordered template of the po- 
larisation modulation produced by the HWP as if it were 
continuously rotated at a mechanical frequency f = ut. 
Continuous rotation of the HWP allows the rejection of all 
the noise components modulated at harmonics different than 
4/ (synchronous demodulation) and is typically employed 
by experim ents optimised to measure the polarisation of the 
CMB (e.g.. I Johnson et al.ll2007l ;l Reichborn-Kiennerud et al.l 
In such experiments, the HWP modulation curve 
leaves a definite synchronous imprint on the time-ordered 
bolometer data streams (timelines), hence it is of utter im- 
portance to characterise the template and remove it from the 
raw data. However, a time-ordered HWP template would be 
of no use to a step-and-integrate experiment such as BLAST- 



Pol, whose timelines are not dominated by the HWP syn- 
chronous signal. 

We have measured the phase shift of the HWP across 
the wavelength range of interest to be ~160°, which appre- 
ciably deviates from the ideal 180°; this is primarily due 
to the unavailability on the market of sapphire substrates 
with the exact desired thickness. However, the modulation 
efficiency of the HWP is only mildly affected by this depar- 
ture from ideality, being above 98% in all three BLASTPol 
bands. Moreover, departures of similar amplitude are not 
uncommon for HWPs at mm and submm wavelengths. 

The goal of the second part of this work was to include 
the measured non-idealities of the HWP as-built in a map- 
making algorithm. We have focused on the implementation 
of a naive binning technique for the case of BLASTPol, un- 
der the assumption of white and uncorrelated noise. As a 
proof of concept, we have presented a preliminary polar- 
isation map for one of the scientific targets observed by 
BLASTPol during its first Antarctic flight, completed in 
January 2011. The inferred direction for the local magnetic 
field in the Carina Nebula star- forming region is in exc ellent 
agreement with the results obtained bv lLi et all l|2006f ) with 
the SPARO instrument. The empirical approach presented 
in this paper will ultimately yield unprecedented accuracy 
on astronomical measurements of the polarisation angle on 
the sky at submm wavelengths. 
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